library(readxl)
install.packages("DESeq2")
options("repos" = c(CRAN="http://mirrors.tuna.tsinghua.edu.cn/CRAN/"))#设置镜像 来于https://zhuanlan.zhihu.com/p/666913516
options(BioC_mirror="http://mirrors.tuna.tsinghua.edu.cn/bioconductor/")
#安装加载BiocManager包
if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)
library(BiocManager)
if (!requireNamespace("BiocManager", quietly = TRUE))#源于https://www.jianshu.com/p/e7c77d16916d
  install.packages("BiocManager")
BiocManager::install("DESeq2")
library(DESeq2)

mRNAdata<-read.csv("E:/R Studio/shixi/shixi 6/mRNA_exprSet.csv",,header = TRUE,sep=",",quote="",row.names=1,fill=TRUE,,check.names=FALSE)#读取csv文件,加入check.names=F，是为了不让列名中的‘-’变成‘.’
metadata<-data.frame(TCGA_id=colnames(mRNAdata)[-1])#把样本名称变成数据框格式,此时数据框中只有TCGA_id 参考https://zhuanlan.zhihu.com/p/554523805
table(substring(metadata$TCGA_id,14,15))#用table函数统计样本的分类.最终结果是有两类：01有多少，11有多少 
sample<-ifelse(substring(metadata$TCGA_id,14,15)=="11","normal","cancer")#TCGA_id第14和15位为11的就是normal，其它就是01异常的cancer
metadata$sample<-as.factor(sample)
View(metadata)#此时metadata中有TAGA_id和sample（即normal和cancer）
countdata<-as.data.frame(mRNAdata[-1])#读取mRNAdata数据去掉首列
colnames(countdata)
library(DESeq2)
dds<-DESeqDataSetFromMatrix(countData=countdata,colData=metadata,design=~sample,tidy=FALSE)#构建数据（表达数据，样本分组信息，模型）
dds2<-DESeq(dds)#差异分析
res<-results(dds2)#得到结果
topGenes<-head(rownames(res[order(res$padj),]),10)#查看前10个按padj值排序的基因的名称
write.csv(res,file="E:/R Studio/shixi/shixi 6/DESeq2_result.csv")#输出差异表达基因
DESeq2_result<-read.csv("E:/R Studio/shixi/shixi 6/DESeq2_result.csv",header = TRUE,sep=",")#读取差异表达基因文件
up<-DESeq2_result[DESeq2_result$pvalue<0.05&DESeq2_result$log2FoldChange>1,]#选择上调的
down<-DESeq2_result[DESeq2_result$pvalue<0.05 & DESeq2_result$log2FoldChange< -1,]#选择下调的
DESeq2_result$Significance<-ifelse(DESeq2_result$pvalue<0.05&DESeq2_result$log2FoldChange>1,"up",
                                   ifelse(DESeq2_result$pvalue<0.05 & DESeq2_result$log2FoldChange< -1,"down","Not_Significant"))
library(ggplot2)
volcano_plot<-ggplot(DESeq2_result,aes(x=log2FoldChange,y=-log10(pvalue),color=Significance))+ #进行相关参数的设置
  geom_point(alpha=0.6,size=2)+
  scale_color_manual(values=c("up"="red","down"="blue","Not_Significance"="grey"))+
  labs(x="log2FoldChange",y="-log10(pvalue)",color="Significance")
volcano_plot #查看图片结果
ggsave("E:/R Studio/shixi/shixi 6/volcano_plot.jpg",volcano_plot)#保存图片

